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Abstract 

We formulate a phase reduction method for a general class of noisy limit cycle oscillators and 
find that the phase equation is parameterized by the ratio between time scales of the noise and 
amplitude-relaxation time of the limit cycle. The equation naturally includes previously proposed 
and mutually exclusive phase equations as special cases. The validity of the theory is numerically 
confirmed. Using the method, we reveal how noise and its correlation time affect limit cycle 
oscillations. 
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Self-sustained oscillations are widely observed in physical, chemical and biological systems 
sl- The oscillations are often described as limit cycle oscillators. Since limit cycle 
oscillators show rich and varied properties, they have been extensively studied as a central 
issue of nonlinear science. Timing of limit cycle oscillation can be described by a single phase 



variable. The phase reduction method is a powerful analytical tool to approximate high- 
dimensional limit-cycle dynamics as a closed equation for only the single phase variable [1 1 . 
Based on the phase description, studies have revealed fascinating properties of limit-cycle 



fla. 



oscillators like response properties and their collective dynamics 

While the theory of phase reduction has been developed mainly for deterministic limit 
cycle oscillators, oscillators in the real world are often exposed to noise. Sources of the 
noise can be internal fluctuations, background noise and also input signals which have noise- 
like statistics f?]. Since noisy limit cycle oscillators also show various nontrivial properties, 
there have been many recent studies of them 0, y, Q, Q, Q, Q, Q, Q. While the 



phase-reduction method is among the most useful ways to study the effects of noise on 
oscillators, two mutually exclusive phase equations have been proposed for a limit cycle 
oscillator driven by white Gaussian noise. The first one is formally the same as the phase 
equation obtained from deterministic oscillators and is in a sense a limiting case of colored 
„o.e wM. ... ... .,,,.0.., ... be.. p.opo.,„„a, to 

square of noise strength and is the technically correct phase equation for white noise jisl . 

Their relationship and which of them is more appropriate description of noisy physical 
oscillators have not been addressed in the literature. Rather, it was recently pointed out 
that both of them fail to describe noisy oscillations in some cases {lo!]. These facts must 
imply existence of a more appropriate phase equation, which will be a starting point for 
future research of noisy oscillations. In this letter, we solve these problems by formulating 
the stochastic phase reduction with careful consideration of relationship between correlation 
time of the noise and relaxation time of the amplitude of the limit cycle. 

n 

Noise in the real world has small but finite correlation time [17[. When the correlation 
time is much smaller than characteristic time scales of the noise-driven system, we can use 
the white noise description by taking the limit where the correlation time goes to zero. For 
limit cycle oscillators, this condition might seem to mean that the correlation time is much 
smaller than the period of oscillation. However, limit cycle oscillators always have other 
significant time scales, i.e., the rate of attraction of perturbations to the limit cycle. These 
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rates characterize stability of the hmit cycle against amplitude perturbation. When the 
limit cycle is very stable to perturbations, the decay time constant could be as small as the 
short correlation time of the noise. Since interplay of small time constants can play a crucial 
role in stochastic dynamical systems, we should carefully consider their relationship when we 
take the white noise limit for noisy limit cycle oscillators. We employ an Ornstein-Uhlenbeck 
process which explicitly has a finite time correlation and then take the white noise limit of 
the process while at the same time keeping track of the time constant for attraction to the 
limit cycle. 

Let us consider a smooth limit cycle oscillator driven by the Ornstein-Uhlenbeck process 
with the time constant r^, 

X = F{X) + aG{X)r]{t) 

(1) 

Tr^V = -r] + ^{t), 

where X{t) G is the state of the oscillator at time t, -F(X) is its intrinsic dynamics, 
G{X) is a vector function, ^(t) is the zero mean white Gaussian noise of unit intensity, 
{^{t)) — ^^"^ {^{^)^{^)) — ^(^ ^ ■5)5 ^^'i then ri{t) represents the zero mean Ornstein- 
Uhlenbeck process with correlation time r,,, {'r]{t)ri{s)) = exp(— |t — s|/r^)/(2r^). As we take 
the limit ^ 0, //(t) approaches the white Gaussian process of unit strength, a represents 
noise strength. F{X) has a stable limit cycle solution Xo(t) satisfying Xq = F{Xq) with 
period T, Xo{t + T) = Xo(t). The phase variable (p is defined around the limit cycle solution 
and increases by T for every cycle of X{t) along the limit cycle. Thus, intrinsic angular 
velocity of the phase is equal to one. We introduce the other — 1 dimensional coordinates 



0- 



(pi,P2, ■ • • ) to describe the dimensional dynamics of X using the coordinate (0, p) 



151]. Without loss of generality, we can shift the origin of p to p = on the limit cycle 
solution. For simplicity of the analysis, we assume that N = 2. Generalization of results 
to any values of A^ is straightforward. We now introduce new variable y(t) = rj{t)^/T^. 



Unlike r/, y has the steady distribution, Po(?/) = exp(— ?/^)/A/7r, which is independent of the 
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correlation time r^. Variable translations from X to {(f), p) and from t] to y gives 

0=1 + ah{(f),p) ^ 



1 y 
P= ^-Tfi4>,p) + (^9{^,P)^ (2) 

The functions h, f and g are defined as h{(f),p) = Vx0 ■ Qi^^)\x=x{(j,,p), f{'t'iP)/^p = 
VxP ■ F{X)\x=x{<t>,p) and 5'(0,p) = VxP ■ G{X)\x=x{<t>,p) \y^- Since the limit cycle at 
p = is stable, we explicitly introduced amplitude-relaxation time of the limit-cycle as Tp, 
which generally depends on (p and assumed that f{(f),0) = and df{(j),0)/dp = —1. The 
value of Tp can be very small if the limit cycle is stiff against amplitude perturbations. 

To eliminate the amplitude variable p and perform the phase reduction, we assume that 
the limit cycle is sufficiently stable and take the limit —>■ 0. Simultaneously, we have to 
take the white noise limit —>■ 0. To consider these two limits at the same time, we take 
the both limits and — simultaneously keeping the ratio k = Tr^/Tp constant. 

Introducing a small parameter e = ^yT^, we translate the variable p to r = p/e, which 
remains 0(1) as e ^ 0. Expanding h, f and g as h{(j), er) = /io(0) + /ii(0)er + /i2(0)e^r^ + . . . , 
/(0,er) = -er + /2(0)eV + /3(0)eV+^^ and^(0,er) = ^o(0) +^i(0)er + ^2(0)eV + . . . , 



we obtain the Fokker-Planck equation 18|, ll9| for the distribution function Q{(f), r, y, t) from 



the stochastic differential equation Eq.(l2]) as 

e'^ = {L,-eL,-e'L,)Q + 0{e'), (3) 

where linear operators are defined as LqQ = {yQ)y + Qyy/'^ + k{rQ)r — crygoQr, LiQ = 
cry[gi{rQ)r + (hoQ)^] + kf2{r^Q)r and L2Q = (ry[g2{r^Q)r + r{hiQ)^] + + kf^{r^Q)r. 
Subscript x means partial derivative with respect to the variable x. We assume that Q 
vanishes rapidly as, y ±00 or r ^ ±00. Expanding Q in a perturbation series in e, 
Q = Qo + eQi + e^<52 + ■ • • , and equating coefficients of equal power of e in Eq. Q, we 
obtain 

e° : L,Q^ = (4) 

: LoQi = LiQo (5) 
d 

: LoQ2 = -^Qo + L2Qo + LiQi. (6) 
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The lowest order equation, Eq. (jlj), has a solution, Qq = P{(f),t)W{(f),r,y), where 
W{(j),r,y) = Vk{l + k) / (agoir) exp{—y'^ — k{y — (1 + k)r / {agojY) is the steady Gaussian 
distribution function of r and y with frozen and g{(j), r) = go{(f)). P{(f), t) is the distribution 
function of the (h. Our primary eoal is to find the evolution equation for P, which is nothiner 
but the reduced Fokker-Planck equation for the phase variable p^, I19I |. 

Since the linear operator Lq has the zero eigenvalue, Eq. ([5]) and ([6]) have to fulfill a 
solvability condition known as the Fredholm alternative. That is, LqU = b has a solution if 
and only if, b is orthogonal to the nuUspace of the adjoint of Lq. This nullspace is simply the 
constant function 1. Thus we can solve LqU = b when the integral of b over {r,y) vanishes. 
To obtain this condition, we integrate both sides of these equations with respect to both 
r and y from — oo to oo. We will see that the condition for Eq. ([6]) is nothing but the 
desired Fokker-Planck equation for (p. Equation ([5]) is solvable since integration over (r, y) 
is zero. To see why, note that integration of the term (rQo)r with respect to r vanishes since 
rQo{r,y) vanishes as |r| — > oo. Integration of yQo{y,r) first with repect to r yields an odd 
function of y which is absolutely integrable and thus its integral over y vanishes. We do 
not need the full expression for Qi at this point, so defer its calculation to the next step. 
Integration of Eq. ([6]) gives 



= Pt 



a 



hoi I {yQi)drdy+ 2{l + k) ^^'^ 



oo J — oo 



+ (7) 



where we used the rapidly vanishing assumption of Q. The coefficient of the 3rd term 
comes from the relationship f^^{yf^W)drdy = ago/ (2 (1 + A;)), which is the correlation 
between y and r for fixed 0. To evaluate J^^{yQi)drdy of the 2nd term, we integrate 
Eq. (j5]) with respect to r from — oo to oo and obtain 

y r Qidr) +UrQidr) =^^!^ye-y\ (8) 

^-OO /y ^ \J~00 J yy 

Since Eq. (|8]) is a differential equation for Q\dr with respect to y, we obtain Q\dr = 
—a {hQP)^ye~y^ / y/n by solving this equation. Then we find that 

/OO /"OO 
/ (yQ,)drdy = --{hoP)^. (9) 
-OO J ~oo ^ 

Substituting Eq. ([9]) into Eq. ([7]) gives the partial differential equation for P as, 

(10) 



= (P* + P^) - - 



{ho{hoP)^)^ - Y^ihigoP)^ 
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which is just the Fokker-Planck equation for the phase variable. Finally, we obtain the phase 
equation as the Ito stochastic differential equation equivalent to the Fokker-Planck equation 

as 

0=1 + + + aZmit), (H) 

where we introduce Z{(f)) = /io(0) = /i(0,O) and = hi{(f))go{(f))/2 = 0)^(0, 0)/2. 

This is also equivalent to the stochastic differential equation 

= ^ + TTkiff'^^'^^ + ^z{mt), (12) 

in the Stratonovich interpretation. 

We now examine the consequence of the above result. The obtained phase equation 
is explicitly parameterized by the ratio between time constants, k = T.rjTp. When the 
correlation time of the noise is much smaller than the decay time constant, we can assume 
A; = and Eq. ( |T2l) is reduced to = l+cr^F(0)+(TZ(0)^(t), which is just the phase equation 
proposed by Yoshimura and Aral . This implies that when noise is white Gaussian noise 
in the strict sense, the 2nd term must be included in the phase equation. On the 

other hand, when the amplitude of the limit cycle decays much faster than the correlation 
time of the noise, or the limit-cycle is sufficiently stable against amplitude perturbations, 
we can assume that k = oo and the 2nd term vanishes. Thus Eq. (\Tp^ is reduced to 



1 + aZ((f))^{t), which is the same to the equation used in j^, 1^, llO, U, Q, Q- The 
latter equation is directly obtained if we apply the standard phase reduction method to 
X = F{X) + aG{X)^{t) without concern for stochastic nature of the perturbation 
Thus, the above result ensures that we can formally use the standard phase reduction in 
these cases. While Eq. f|T2l) agrees with previously proposed equations at opposite limits of 
the parameter k, it deviates from both of them in the middle range of k. Therefore, we can 
conclude that in order to properly describe stochastic phase dynamics for a general value 
of k, we must consider the coefficient of the 2nd term correctly as + k) in the phase 
equation. 

To see the effect of the weight 1/(1 + we will calculate the steady distribution function 
for the phase. Requiring the steady condition = to Eq. (fTOj) . we obtain the steady 
distribution as: 



Po(0) = ^(l+^' 



2 TTW) ' 



0{a% (13) 
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where we used power series expansion of the distribution in terms of o"^. Qq is defined as 
^0 = T^^ J^^ Y {(j)) / {1 + k{(j)))d(j). As we increase noise strength a from zero, the phase 
distribution starts to deviate from 1/T of non-perturbed oscillators. While magnitude of 
the deviaton is a function of a, actual shape of this depends on the ratio k{(j)). 

Using the steady distribution, we can calculate the mean frequency of the noisy oscillator 
defined as = \imt-,oct~^ Jq <P{t)dt. Replacing the long term average with the ensemble 



average, i.e. 1) = <j)PQ{(f))d(f), and substituting the Ito equation Eq. ffTTl) into <j), we have 

n = l + a^no + 0{a^), (14) 
where we used the fact that (f){t) is independent from ^(t) in the Ito equation. As pointed out 



in the previous study [15|, the mean frequency depends on the noise strength. In addition 
to the strength, our result reveals that the frequency also depends on and Tp through the 
ratio k. As we change these values, the mean frequency will increases or decreases depending 
on the sign of Qq- 

In order to validate the above analysis, we numerically examine stochastic phase dy- 
namics and calculate Pq and Q directly from the stochastic differential equation ([1]). 
As a simple example, we use the Stuart-Landau (SL) oscillator, X = {x,y), F{X) = 
{^{Z{W)),'^{Z{W))), where W = x + iy and Z{W) = {X{l+ic) +iuj)W - X{l + ic) \W\^W, 
which is rescaled such that amplitude relaxation time will explicitly appear. We define 
phase and amplitude coordinates (0, r) as = (arctan(y/x) — clog(x^ + y'^)/2)/uj and 
r = a/x^ + — 1. The limit cycle solution + = 1 is given as r = in the coor- 
dinate. The decay time constant to the limit cycle solution is Tp = 1/(2A). Figure 1 shows 
steady state distributions of the phase for various values of time constants and Tp. As 
expected, the distribution changes as a function of time constants. Distributions, however, 
are the same as far as the ratio between them is the same. Numerical results are well fitted 
by the analytical result Eq. (fT3|) . Figure 2 shows the mean frequency Q a.s a function of r,, 
and Tp. As indicated by the above analysis, Q increases as a function of and decreases 
as a function of Tp. Theoretical predictions, Eq. (fT4l) . agree fairly well with the numerical 
results. 

The above results clearly indicate that, when we eliminate fast variables in stochastic 
dynamical systems, characteristic time scales of the fast variables should be seriously con- 
sidered even though variables themselves are eventually eliminated. In particular, white 
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Gaussian noise is actually an idealization of physical processes with small but finite time 
correlation. Interactions between small time scales can give crucial effects to stochastic dy- 
namics. Thus similar situations may also arise even when we use reduction methods other 
than the phase reduction to stochastic phenomena j2o|. Actually a similar situation arises 



in the analysis of classical Brownian motion with inertia 21|. The above results also tell us 
that dynamical systems driven by the white- Gaussian noise are derived through reduction 
methods not only from literally white-noise-driven systems but also from systems driven 
by realistic noise with finite time correlations. The non-agreement between previously pro- 
posed phase equations is due to this ambiguity. Our results ensure that we can choose the 
most suitable reduced equation as far as we explicitly indicate time scales of the noise and 
dynamical systems. 

In summary, we have formulated stochastic phase reduction for a general class of smooth 
limit cycle oscillators. The derived stochastic phase equation is parameterized by the ratio 
between the correlation time of the noise and the decay time of amplitude perturbations. 
Whereas previously proposed phase equations are realized only at opposite limits of the 
ratio, the obtained phase equation is valid in the whole range of values of the ratio. We 
have calculated steady phase distributions and the mean frequency of the noisy oscillator 
and reveal their dependence on the time scales. The results suggest significance of fast time 
scales in reduction methods of stochastic phenomena. 
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